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Glycogen synthase kinase-3 (Gsk-3) activity is an important regulator of numerous signal 
transduction pathways. Gsk-3 activity is the sum of two largely redundant proteins, Gsk- 
3a and Gsk-3fi, and in general, Gsk-3 is a negative regulator of cellular signaling. Genetic 
deletion of both Gsk-3a and Gsk-3fi in mouse embryonic stem cells (ESCs) has previously 
been shown to lead to the constitutive activation of the Wnt/p-catenin signaling pathway. 
However, in addition to Wnt signaling, all Gsk-3-regulated pathways, such as insulin sig- 
naling, are also affected simultaneously in Gsk-3ar^^; ESCs. In an effort to 
better understand how specific signaling pathways contnbute to the global pattern of gene 
expression in Gsk-3a^^^; Gsk-3^^^^ ESCs, we compared the gene expression profiles in 
Gsk-3a^^^; Gsk-3^^^^ ESCs to mouse ESCs in which either Wnt/p-catenin signaling or 
phosphatidylinositol 3-kinase (PI3K)-dependent insulin signaling are constitutively active. 
Our results show that Wnt signaling has a greater effect on up-regulated genes in the Gsk- 
3ar^^] Gsk-3^^^^ ESCs, whereas PI3K-dependent insulin signaling is more responsible 
for the down-regulation of genes in the same cells. These data show the importance of 
Gsk-3 activity on gene expression in mouse ESCs, and that these effects are due to the 
combined effects of multiple signaling pathways. 

Keywords: embryonic stem cells, glycogen synthase kinase-3, Wnt, phosphatidylinositol 3-kinase, microarray, 
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INTRODUCTION 

Glycogen synthase kinase-3 (Gsk-3) is an intracellular ser- 
ine/threonine kinase activity that was originally identified as the 
rate-limiting step in the glycogen synthesis pathway (1). When the 
genes encoding Gsk-3 were cloned, Gsk-3 activity was revealed 
as the combined action of two distinct genes, Gsk-3a and Gsk- 
3P (2). Gsk-3 isoforms show a high degree of similarity at the 
amino acid level, including 98% similarity within the catalytic 
domain. Despite divergence at the amino- and carboxyl-termini 
of each isoform, it appears that the unstructured carboxyl-termini 
are essential for Gsk-3 activity (3). Rare among kinases, Gsk-3 is 
active at a basal state, while pathway activation from upstream 
signaling cascades results in the inhibition of Gsk-3 activity (4-6). 

There are two primary mechanisms responsible for Gsk-3 inhi- 
bition - one that regulates the phosphorylation of Gsk-3a and 
Gsk-3p on serines 21 and 9, respectively, and another that is 
independent of this phosphorylation event. Several signaling path- 
ways, including protein kinase A (PKA), Hedgehog, transforming 
growth factor-p (TGF-fi), nuclear factor of activated T-cells (NF- 
AT), and phosphatidylinositol 3-kinase (PI3K)-dependent insulin 
signaling, all affect Gsk-3 activity via phosphorylation of the 



amino-terminal serines (4-6). Insulin binding to its receptor trig- 
gers the activation of PI3K, which phosphorylates and activates 
Akt, which in turn phosphorylates and inhibits Gsk-3 activity (7). 
We have previously shown that stable expression of a constitutively 
active form of the pi 10a catalytic subunit of PI3K (termed pi lO*) 
(8) in WT mouse ESCs can effectively lead to the phosphorylation 
of Gsk-3 (9), and expression of pi 10* has been shown to activate 
insulin signaling in mouse ESCs (10). 

Gsk-3 activity also plays a central role in Wnt signaling as it is 
part of a cytosolic protein complex termed the f5-catenin destruc- 
tion complex (11). The role of Gsk-3 is to directly phosphorylate 
P-catenin (12,13), which targets P-catenin for ubiquitin-mediated 
degradation ( 14), keeping the Wnt pathway inactive. Activation of 
the Wnt signaling pathway inhibits Gsk-3 activity, but this effect 
is not mediated through amino-terminal serine phosphorylation. 
Instead, upon Wnt ligand binding to a Frizzled/Lrp receptor com- 
plex, the destruction complex is neutralized via the translocation 
of Gsk-3 to the cell membrane, where it phosphorylates Lrp (15, 
16). Gsk-3 is then subsequently endocytosed into multivesicular 
bodies, further insulating it from P-catenin (17). Once P-catenin 
protein levels accumulate to high levels in the cytoplasm, it is 
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subsequently translocated to the nucleus where it can directly com- 
plex with lymphocyte enhancer factor (Lef)/T-cell factor (Tcf) 
proteins to activate Wnt target genes (18, 19). Gsk-3a and Gsk- 
3P are redundant with respect to their roles in regulating Wnt 
signaling; the genetic deletion of either Gsk-3a or Gsk-3^ alone 
is insufficient to activate Wnt signaling (20-25). Only upon the 
deletion of three of the four total alleles for Gsk-3a and Gsk-3fi is 
Wnt pathway activation achieved (23). Double knockout of both 
Gsk-3a and Gsk-3p in mouse embryonic stem cells (ESCs) (23) 
or mouse neural progenitors (26) results in the constitutive acti- 
vation of the Wnt pathway. Importantly, while Gsk-3 regulates 
both insulin and Wnt signaling, these pathways do not activate 
one another; Wnt signaling does not activate insulin signaling and 
insulin signaling does not activate Wnt signaling (27-29). 

Lithium is a direct inhibitor of Gsk-3 activity (30, 31) and is 
used therapeutically for the treatment of bipolar disorder (32). The 
inhibition of Gsk-3 activity has also been shown to be a potential 
means to reduce the pathology associated with Alzheimer's disease 
(33-38). Because the dysregulation of Gsk-3 activity is believed 
to be an important contributor to many diseases, including 
Alzheimer's disease, cancer, mental illness, and diabetes, the acqui- 
sition of more fundamental knowledge about the effects of reduced 
Gsk-3 activity could be important in the context of these diseases. 

In this paper, we describe the analysis of genome-wide gene 
expression patterns in Gsk-3a^^^; Gsk-3^^^^ ESCs. Since the loss 
of Gsk-3 activity results in the activation of numerous signaling 
pathways simultaneously, we further dissected the effects of los- 
ing Gsk-3 activity by examining gene expression patterns in ESCs 
with either constitutively active PI3K-mediated insulin signaling 
or constitutively active Wnt signaling. 

MATERIALS AND METHODS 
MOUSE EMBRYONIC STEM CELL CULTURE 

wild-type (WT) mouse ESCs (E14K), Gsfc-3a"^"; Gsfc-jp^^" (23), 
pi 10'' (9), and P-catenin S33A (39) stable cell lines, along with 
respective control WT ESCs stably transfected with puromycin- 
resistant plasmids, were grown on 0.1% gelatin-coated plates in 
high glucose DMEM (Life Technologies) supplemented with 15% 
fetal bovine serum (Hyclone), 1% non-essential amino acids, 1% 
sodium pyruvate, 1% L-glutamine, 1% penicillin/streptomycin 
(Life Technologies), 55 p.M 2-mercaptomethanol, and 1000 U/ml 
ESGRO (Millipore). Media was replenished every other day. 

RNA ISOLATION, cDNA SYNTHESIS, AND QUANTITATIVE PCR 

RNA was isolated from 5 x 10^ to 1 x 10*^ ESCs using the 
MirVana Total RNA Isolation Kit (Applied Biosystems) accord- 
ing to the manufacturer's instructions. Two micrograms of 
RNA were then used for cDNA synthesis using the high 
capacity cDNA reverse transcription kit (Applied Biosystems) 
following the manufacturer's protocol. Quantitative RT-PCR 
was performed on an Applied Biosystems StepOne machine 
using TaqMan master mix and one of the following Taq- 
Man assays (Apphed Biosystems): Axin2 (Mm00443610_ml), 
Brachyury (Mm00436877_ml), Bhmtl (Mm04210521_gl), 
Bhmt2 (Mm00517726_ml), Cdx2 (Mm01212280_ml), Ido2 
(Mm00524206_ml),Anxa8 (Mm00507926_ml), Aqp8 (Mm0043 
1846_ml),Gata6(Mm00802636_ml),Wnt6 (Mm00437351_ml), 



Acta2 (Mm01546133_ml), or Foxql (Mm01157333_sl). Three 
biological replicates and three technical replicates were used for 
each cell type. All threshold cycle (Ct) values were normalized to 
a mouse Gapdh endogenous control (Mm99999915_gl) (Applied 
Biosystems), and relative quantification was calculated from the 
median Ct value (40). P-values were calculated using Data Assist 
software (Applied Biosystems). 

MICROARRAY ANALYSIS 

RNA was isolated using the MirVana Total RNA Isolation Kit 
(Applied Biosystems) according to the manufacturer's instruc- 
tions. The concentration of RNA was measured in a Nanodrop 
ND-1000 UV-Vis spectrophotometer (Nanodrop) and an Agilent 
2100 Bioanalyzer Lab-On-A-Chip Agilent 6000 Series II chip (Agi- 
lent) was used to determine the integrity of the samples. RNA 
from each cell line was hybridized onto a SurePrint G3 Mouse 
GE 8 X 60 K Microarray, 8 x 60 K, AMADID 028005 (Agilent). 
Hybridization was performed overnight at 45°C. We performed 
arrays for each cell line using RNA that was isolated in biological 
triplicate (n = 3). For Gsk-3 DKO ESCs, we used « = 4 biological 
replicates. SurePrint arrays were scanned with an Agilent G2505C 
Microarray Scanner (Agilent). The information about each probe 
on the arrays was extracted from the image data using Agilent 
Feature Extraction 10.9 (EE) and .txt files were generated. The raw 
intensity values from these files in imported into the mathematical 
software package "R," which was used for aU data input, diagnos- 
tic plots, normalization, and quality checking steps of the analysis 
process using scripts developed by Dr. Peter White, Director of 
the Biomedical Genomics Core at Nationwide Children's Hospital 
specifically for this analysis. These scripts call on several Biocon- 
ductor packages, an open-source and open-development software 
project to provide tools for the analysis and comprehension of 
genomic data (41). The algorithm used for normalization of gene 
expression was designed for use with Agilent One-Color Analy- 
sis. The median green (Cy3) intensities are normalized between 
the arrays using the Quantile Normalization package in "R" (42). 
Quantile normalization is a non-linear probe-level normalization 
that results in the same empirical distribution of intensities for 
each array. This is a significantly more robust approach than sim- 
ply normalizing to the median value of each array. The genes that 
were altered by twofold either way and had a false discovery rate 
(FDR) of <10% were sorted and used for further interpretation 
of the microarray data. 

STATISTICAL ANALYSIS OF MICROARRAY DATA 

Statistical analyses were performed using two well-validated and 
commonly used approaches - significance analysis of microar- 
rays (SAM) and adjusted p-value. SAM is a powerful tool for 
analyzing microarray gene expression data useful for identifying 
differentially expressed genes between two conditions (43). SAM 
calculates a test statistic for relative difference in gene 5 expression 
based on permutation analysis of expression data and calculates a 
FDR using the q-value method presented in Ref (44). In outline, 
SAM identifies statistically significant genes by carrying out gene 
specific (--tests and computes a statistic for each gene, which mea- 
sures the strength of the relationship between gene expression and 
a response variable. This analysis uses non-parametric statistics. 
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since the data may not follow a normal distribution. The response 
variable describes and groups the data based on experimental con- 
ditions. In this method, repeated permutations of the data are used 
to determine if the expression of any gene is significant related to 
the response. The use of permutation-based analysis accounts for 
correlations in genes and avoids parametric assumptions about 
the distribution of individual genes. For this experiment, SAM 
analysis was implemented in R using the Bioconductor Siggenes 
package. In a one-color experimental design, a two-class unpaired 
analysis is typically performed for each experimental comparison, 
whereas in a two-color approach a one-class analysis is used. Typi- 
cally, an FDR cutoff in the range of 10-20% is chosen to maximize 
sensitivity without significantly impacting accuracy. For the cur- 
rent study, a 10% FDR was used to generate the list of significantly 
differentially expressed genes. 

Adjusted p-value 

The moderated t-statistic (f) is computer for each probe and for 
each contrast in the experimental design. This has the same inter- 
pretation as an ordinary f-statistic except that the standard errors 
have been moderated across genes, i.e., shrunk toward a com- 
mon value. This has the effect of borrowing information from the 
ensemble of genes to aid with inference about each individual gene. 
The p-value is obtained from the distribution of the moderated 
f-statistic. Finally these p-values are adjusted for multiple test- 
ing using the Benjamini and Hochberg's (45) step-up method for 
controlling the FDR. This is the most popular method for p-value 
adjustment. If all genes with p-value below a threshold, say 0.05, 
are selected as differentially expressed, then the expected propor- 
tion of false discoveries in the selected group is controlled to be less 
than the threshold value, in this case 5%. For the current study, the 
adjusted p-values were calculated using the Bioconductor limma 
package. 

GENE ENRICHMENT ANALYSIS 

After filtering gene lists to remove duplicates, the HGNC symbols 
for genes increased or decreased by twofold or more were entered 
into ToppFun, within the ToppGene Suite (toppgene.cchmc.org) 
(46). ToppGene then performs functional enrichments, looking 
for sets of co-regulated genes. Calculations are performed using 
FDR for correction, and a p-value cutoff of <0.05. 

RESULTS 

We had previously analyzed the genome-wide expression in Gsk- 
3a^^^; Gsfc-Jp^^^ ESCs (hereafter referred to as Gsk-3 dou- 
ble knockout; Gsk-3 DKO) using Affymetrix microarrays, and 
found hundreds of genes whose expression was significantly 
increased or decreased compared to WT ESCs (9). Here, we 
repeated this microarray analysis using the Agilent platform. 
For all of the microarray gene expression analyses, we used a 
twofold or greater change in gene expression as our thresh- 
old for significance. We also isolated RNA and performed the 
hybridization to the microarrays at the same time to reduce 
experiment-to-experiment variability. Consistent with our pre- 
vious results using the Affymetrix platform, we found 1313 genes 
up-regulated twofold or more in Gsk-3 DKO ESCs, while 2178 
genes were down-regulated twofold or more (a complete list of 



genes can be found in Datasheet 1 in Supplementary Material). 
One of the important signaling pathways regulated by Gsk-3 is 
the Wnt pathway (47). We therefore expected our microarray 
data from Gsk-3 DKO ESCs to reveal many of the same Wnt 
target genes that have been demonstrated experimentally in a vari- 
ety of model systems [Wnt Homepage, www.stanford.edu/group/ 
nusselab/cgi-bin/wnt/target_genes]. A few well-established direct 
targets of Wnt signaling, such as Brachyury (T) and Axin2, were 
found to be increased substantially in GsA:-3 DKO ESCs (162.5- 
and 8.6-fold, respectively) (Table 1). We were surprised, however, 
to find that Brachyury and Axinl were the exception of the 65 Wnt 
putative target genes; we identified in our microarray data from 
Gsk-3 DKO ESCs that only seven additional genes (Sp5, Cdxl, 
Stra6, Lefl, Cyclin Dl, PTTG, and Fgfl8) had their expression 
increased more than twofold (Table 1 ) . In fact, 1 3 Wnt target genes 
showed decreased expression of more than twofold (Table 1). 

The relative paucity of up-regulated Wnt target genes was sur- 
prising, especially, since it has been shown that a reporter construct 
containing multimerized Lef/Tcf binding sites is strongly activated 
in Gsk-3 DKO ESCs [(23); unpublished observation]. Therefore, 
we wanted to compare the gene expression data from Gsk-3 DKO 
ESCs with another mouse embryonic stem cell line that has consti- 
tutively active Wnt signaling. WT ESCs stably expressing a form of 
fi-catenin in which serine 33 has been mutated to an alanine (S33A) 
that prevents the phosphorylation by Gsk-3 on this residue as well 
as subsequent ubiquitination, resulting in the constitutive activa- 
tion of Wnt signaling. These cells have previously been shown to 
potently activate the expression of several Wnt target genes (39). 
The P-catenin S33A cells and their control cells were included in 
our microarray experiment to investigate how similar the patterns 
of gene expression were compared to Gsk-3 DKO ESCs. 

Microarray data showed that 1468 genes were up-regulated 
twofold or more in S33A cells compared to control cells, while 
1412 genes were down-regulated twofold or more. As expected, we 
confirmed the high levels of Brachyury, Axin2, and Cdxl that were 
also observed in the Gsk-3 DKO ESCs; however, only three addi- 
tional genes (Sp5, Cdx4, and VEGF) showed a twofold or greater 
increase in gene expression in the S33 A ESCs (Table 1 ) . And similar 
to the data from GsA:-3 DKO ESCs, 12 Wnt target genes are down- 
regulated twofold or more in S33A cells (Table 1). Taken together, 
these data show that some Wnt target genes, such as Brachyury 
and Axinl, are indeed transcribed at high levels in two different 
cell lines in which Wnt signaling is constitutively active - Gsk-3 
DKO and fi-catenin S33A. Our data also show that many Wnt 
target genes are not activated in these cell lines as was predicted. 

Initially, we had expected that many of the genes whose expres- 
sion was significantly increased in Gsk-3 DKO ESCs would be due 
to activation of the Wnt pathway. Since we did not observe the 
increased expression of Wnt target genes, we asked whether the 
activation of another Gsk-3-dependent signaling pathway could 
explain the data. Because the genetic deletion of Gsk-3 isoforms 
likely leads to the activation of signaling pathways other than the 
Wnt pathway, we speculated that the combined effect of multiple 
pathway activation likely leads to new patterns of gene expression. 
Another important signal transduction pathway that is regulated 
by Gsk-3 activity is the insulin signaling pathway. To investi- 
gate the contribution of insulin signaling to the effects on gene 
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Table 1 | Comparison of known Wnt target genes to their expression 
levels in Gsk-3 DKO and S33A ESCs as determined by microarray. 

Gene DKO S33A 



(Continued) 



Table 1 | Continued 



Gene DKO S33A 



MMP2, MMP9 -1.9 1.6 

LBH -2 1.4 

Pitx2 -2.1 -1.5 

ITF-2 -2.3 -1.9 

Twist -2.4 -2 

Periostin -2.4 1 

VEGF -2.5 2.4 

Snail -2.8 -1.5 

Tcf-1 (Hnfia) -3 -1.1 

n-myc -3.4 -3.4 

Gremlin —3.5 —3.5 

EOF receptor -3.8 -1.5 

Delta-like 1 -9.8 -4.3 

WISP -10 -1.1 



The list of genes was found at The Wnt Homepage website (www.stanford.edu/ 
group/ nusselab/ cgi-bin/ wnt/ target_genes). Genes whose expression increased 
more than twofold are shown in red, while genes whose expression decreased 
more than twofold are shown in green. 

expression seen in the Gsk-3 DKO ESCs, we analyzed genome-wide 
microarray gene expression data from WT ESCs stably expressing 
a constitutively active form of the pllO subunit of PI3K (pllO*) 
(9, 10). Eight hundred sixty-four genes were up-regulated twofold 
or more in pi 10* ESCs compared to control cells, while 1660 genes 
were down-regulated twofold or more. 

After compiling the list of genes whose expression was signif- 
icantly changed in Gsk-3 DKO, pi 10*, and S33A ESCs, we were 
curious as to the similarities and disparities in the patterns of 
gene expression in these cell lines. To perform this evaluation, 
we removed all duplicate probes from our lists of genes from the 
microarray experiment, and we then sorted based on changes seen 
in each list. The Gsk-3 DKO and pllO* ESCs had 1313 and 864 
genes, respectively, that were up-regulated twofold or more. Of 
these, 206 genes were up-regulated in both cell types (Figure lA). 
Similarly, 1468 genes were up-regulated by at least twofold in S33A 
cells, and 336 of these genes overlapped with the 1313 genes up- 
regulated in the Gsk-3 DKO ESCs (Figure lA). As we expected, 
there was relatively little overlap in the genes up-regulated in both 
pllO* and S33A ESCs; only 119 genes were in common between 
these gene sets (Figure lA). This comparison of up-regulated 
genes shows a greater degree of overlap between Gsk-3 DKO and 
S33A ESCs than with Gsk-3 DKO and pi 10* ESCs. 

Interestingly, the similarities between gene sets were even 
greater when examining the genes whose expression was down- 
regulated twofold or more. Two thousand one hundred seventy- 
eight genes showed decreased expression in Gsk-3 DKO ESCs, 
compared with 1660 genes in pi 10* ESCs and 1412 genes in 
S33A ESCs. Furthermore, we found that 987 genes were in com- 
mon between pi 10* and Gsk-3 DKO ESCs, while only 516 genes 
were shared between S33A and Gsk-3 DKO ESCs (Figure IB). 
This pattern is striking in that it is the opposite of what was seen 
with the up-regulated genes, with a greater overlap between Gsk-3 
DKO and pi 10* ESCs. In addition, the overlap between pi 10* and 



Brachyury 162.5 41.7 

SP5 75.8 15.1 

Cdxl 16.5 82.5 
Axin-2 8.6 4.6 
Stra6 4.8 1.3 

LEF1 4 -7 

Cyclin D 3.1 -3.5 

Pituitary tumor transforming gene (PTTG) 2.2 1.8 

FGF18 2.1 -2.9 

BMP4 1.8 1.2 

PPAR8 1.6 1.1 

Endothelin-1 1.5 1.4 

Telomerase 1.5 1.4 

LGR5/GPR49 1.5 -1.2 

Jagged 1 .4 1 

sFRP-2 1.4 -1.5 

CD44 1.3 1.4 

c-myc binding protein 1.3 -1.9 

Msll 1.3 1.2 

Nitric oxide synthase 2 1 .3 1 .4 

Nkx2.2 1.3 1.2 

c-myc 1.2 -3.3 

Met 1.2 -1.2 

Claudin-1 1.1 1 

Id2 1.1 1.2 

FoxNI 1.1 1.1 

Gbx2 1.1 -2.2 

MMP-7 1 1 

Osteoprotegerin 1 1 

Wnt3a 1 1 .2 

Neurogenin 1 1 1.1 

Tiami -1.1 -1.2 

FGF20 -1.1 -1.4 

Cdx4 -1.1 2.3 

Sox17 -1.2 -1.6 

Oct-4 -1.2 -1.3 

Runx2 -1.3 1 

Sox2 -1.3 1.4 

Frizzled 7 -1.3 -2.1 

Gastrin -1.4 1.3 

Sox9 -1.4 -1.1 

CCN1/Cyr61 -1.4 -1.2 

Follistatin -1.4 1.2 

NeuroDI -1.5 -5.6 

c-jun -1.6 -2.2 

Nr-CAM -1.6 -1.5 

SALL4 -1.6 -2.3 

Nanog -1.6 1.4 

FGF9 -1.7 1 

lrx3andSix3 -1.8 -1.4 

Cacnaly -1.8 -1.8 
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FIGURE 1 I Venn diagram representing the overlap in genes 
whose expression is changed at least twofold in Gsk-3 DKO, 
pllO*, and fl-catenin S33A ESCs. (A) Representation of genes 
wliose expression was increased in Gsk-3 DKO (blue; 1313 genes 
total), pllO* (red; 864 genes total), and S33A ESCs (green; 1468 
genes total). Genes whose expression is increased in multiple genes 



sets are denoted by the numbers at the interface of the gene sets. 
(B) Representation of genes whose expression was decreased in 
Gsk-3 DKO (blue; 2178 genes total), pllO* (red; 1660 genes total), 
and S33A ESCs (green; 1412 genes total). Genes whose expression 
is decreased in multiple genes sets are denoted by the numbers at 
the interface of the gene sets. 



S33A ESCs was also more pronounced, with the same 289 genes 
down-regulated at least twofold in both cell lines (Figure IB). 
These data suggest that the loss of Gsk-3 in ESCs results in the up- 
regulation of mostly distinct genes in either pi 10* or S33A ESCs, 
but that those genes that are down-regulated in the absence of 
Gsk-3 share a greater overlap in cells with activated PI3K signaling 
and activated Wnt signaling. 

Since different sets of genes were either increased or decreased 
in expression in pllO* and S33A ESCs, we bioinformatically per- 
formed enrichment analyses to see if functionally related genes 
were co-regulated in a signaling pathway-specific fashion (Tables 2 
and 3). Genes that were up-regulated in Gsk-3 DKO ESCs were 
involved in gene ontology (GO): molecular function of protein 
domain (GO:0019904) and also the GO: biological process for cil- 
ium organization (GO:0044782). In addition, genes encoding for 
CD molecules, tumor necrosis factor receptor superfamily, and 
intraflagellar transport homologs were up-regulated. Genes that 
were increased in expression in p 1 1 0* ESCs showed an enrichment 
for the GO: molecular function of piRNA binding (GO:0034584) 
and an showed an increase in expression in genes encoding for zinc 
fingers, C2H2, and BTB domain-containing (ZBTB). In addition, 
both Gsk-3 DKO and pi 10* ESCs showed increased expression 
of genes involved in urothelium and lower urinary tract devel- 
opment. For genes up-regulated in S33A ESCs, the top GO: 
molecular function clusters were genes involved in organic acid 
binding (GO:0043177) and protein homodimerization activity 
(GO:0042803). Also, a significant number of genes that contain 
homeobox domains were up-regulated in the S33A ESCs, as well 
as a high number of genes encoding for CD molecules. 

We found it notable that the enrichments among down- 
regulated genes were more robust than among the up-regulated 
genes (Table 4). For genes decreased twofold or more in Gsk-3 
DKO ESCs, the greatest enrichment for GO: molecular function 
was in genes encoding for sequence-specific DNA binding, i.e., 
transcription factors (GO:0043565). In addition, a large number 



of genes involved in the development of the cardiovascular sys- 
tem (GO:0072358) were down-regulated in Gsk-3 DKO ESCs. 
Also decreased in expression were genes involved in extracellu- 
lar matrix organization and elastic fiber formation. For down- 
regulated genes in pi 10* ESCs, the most enriched GO: molecular 
function was for genes involved in receptor binding (GO:0005 102), 
while the highest enrichment for GO: biological processes were for 
organ morphogenesis (GO: 0009887), extracellular matrix orga- 
nization (GO:0030198), and cardiovascular system development 
(GO:0072358). Furthermore, many genes involved in extracel- 
lular matrix organization were decreased in pi 10* ESCs. Genes 
down-regulated in S33A ESCs were enriched for genes encoding 
sequence-specific DNA binding (GO:0043565), as well as enrich- 
ment for genes involved in neurogenesis (GO:0022008) and car- 
diovascular system development (GO:0072358). The results from 
these enrichment analyses showed several similarities and differ- 
ences with respect to the genes that have increased or decreased 
expression in Gsk-3 DKO, pi 10*, and S33A ESCs, and provide 
a framework for beginning to better understand the complex 
interrelationships between PI3K-mediated signaling and Wnt sig- 
naling, as well as providing important insights into the patterns of 
gene expression seen in Gsk-3 DKO ESCs. 

While the data obtained from microarray experiments can 
provide an excellent overall view of differential patterns of gene 
expression, it is nonetheless very important to use an independent 
technique to experimentally validate the changes in gene expres- 
sion that were observed via microarray. Therefore, we selected a 
handful of genes whose expression patterns we of interest, and we 
performed quantitative reverse-transcriptase PGR (qPGR) using 
TaqMan probes on RNA isolated from Gsk-3 DKO, pi 10*, S33A 
ESCs, and their respective control cell lines. As an initial quality 
control measure, we performed qPCR to validate the expression of 
the known Wnt target genes, Axinl and Brachyury. As expected, 
we saw substantial increases in gene expression in both Gsk-3 DKO 
and S33A ESCs for Axinl (up 7-fold in Gsk-3 DKO and 3.7-fold 
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Table 2 | Summary of enrichment analysis of genes up-regulated in 
Gsk-3 DKO, pllO*, and S33A ESCs. 



Name ID p-Value Hit count 



Table 3 | Summary of enrichment analysis of genes down-regulated in 
Gsk-3 DKO, p110*, and S33A ESCs. 



Name ID p-Value Hit count 



GO: molecular function 



Protein domain specific binding 


GO:0019904 


1.01 X 10" 


-5 


55/633 


GO: biological process 










Cilium organization 


GO:0044782 


6.06X 10" 


-6 


19/123 


Gene family 










CD molecules 


CD 


5.14 X 10" 


-8 


21/276 


Tumor necrosis factor receptor 


TNFRSF 


5.77 X 10" 


-5 


4/12 


superfamily 










Intraflagellar transport homologs 


IFT 


1.52 X 10" 


-4 


4/15 


Coexpression Atlas 










Mendel_RNAseq_e17.5_urothelium 




2.74x 10- 


16 


101/989 


Developing lower urinary 




1.04x 10" 


13 


82/802 


tract_el3.5 










P110'' ESCs ^^^^^^^^^^^^^^^^^^^^^^^H 


GO: molecular function 










piRNA binding 


GO:0034584 


3.17 X 10" 


-5 


3/3 


Gene family 










Zinc fingers, C2H2, and BTB 


ZBTB 


5.33 X 10- 


-5 


6/47 


domain containing 










Coexpression Atlas 










Mendel_RNAseq_e17.5_urothelium 




1.58 X 10- 


10 


65/989 


Developing lower urinary 




1.30 X 10" 


-6 


30/802 


tract_e13.5 










GO: molecular function 










Organic acid binding 


GO:0043177 


4.93x lO- 


-6 


26/230 


Protein homodimerization 


GO:0042803 


S.ISx 10" 


-5 


56/674 


activity 










Gene family 










CD molecules 


CD 


5.08 X 10- 


-5 


18/276 



Bonferroni correction method. 
Al! p-Value cutoffs - 0.05. 



in S33A ESCs) and Brachyury (up 622-fold in Gsk-3 DKO and 
38-fold in S33A ESCs), while both genes has reduced expression 
in pi 10* ESCs (0.4- and 0.7-fold changes, respectively) (Figure 2). 
These results are consistent with Axin2 and Brachyury being Wnt 

target genes. 

We then proceeded to validating the expression of addi- 
tional genes whose expression levels were shown to be substan- 
tially increased or reduced in the microarray experiments. We 
selected four genes, Bhmtl (Betaine-homocysteine methyltrans- 
ferase), Bhmt2 (Betaine-homocysteine methyltransferase 2), Cdx2 
(Caudal-type homeobox2), and Ido2 (Indoleamine 2,3-dioxygenase 
2) , whose expression, by microarray, was shown to be up-regulated 



GO: molecular function 



Sequence-specific DNA binding 


GO:0043565 


8.59 X 10- 


-13 


109/726 


Sequence-specific DNA binding 


GO:0003700 


728 X 10" 


-12 


141/1067 


transcription factor activity 










GO: biological process 










Cardiovascular system 


GO:00722358 


2.15X 10" 


-21 


148/871 


development 










to 11 1 wciy 










Extracellular matrix organization 




1.40 X 10- 


-9 


50/264 


Elastic fiber formation 




1.51 X 10" 


-7 


15/41 


GO: molecular function 










Receptor binding 


GO:0005102 


1.51 X 10" 


-11 


133/1341 


GO: biological process 










WigdM 1 1 lUI pi luyciicoib 


O*— '.UUUooo / 


9 nq w -l^- 


-34 


1 /L7 /07R 
1 / /O /o 


Extracellular matrix organization 


GO:0030198 


1.59x 10" 


-22 


70/330 


Cardiovascular system 


GO:0072358 


2.34x 10" 


-22 


124/871 


development 










Pathway 










Extracellular matrix organization 




6.17 X 10" 


-16 


54/264 


GO: molecular function 










Sequence-specific DNA binding 


GO;0043565 


2.15x 10" 


-10 


81/726 


Sequence-specific DNA binding 


GO:0003700 


3.54 X 10" 


-10 


105/1067 


transcription factor activity 










GO: biological process 










Neurogenesis 


GO:0022008 


3.88 X 10" 


-15 


143/1360 


Extracellular matrix organization 


GO:0072358 


1.21 X 10" 


-14 


104/871 



Bonferroni correction method. 
All p-Value cutoffs - 0.05. 



in both Gsk-3 DKO and S33A ESCs, but not p 1 1 0* ESCs. Validation 
by qPCR showed that the expression of Bhmtl was up 14.1- and 
70.1-fold in Gsk-3 DKO and S33A ESCs, respectively, compared 
to a modest 2.8-fold increase in pi 10* ESCs (Figure 3). Similarly 
Bhmt2 expression was shown to be increased 13.1-fold in Gsk-3 
DKO ESCs and 56.5-fold in S33A ESCs, whUe expression was only 
increased by 1.3-fold in pi 10* ESCs (Figure 3). The expression 
of Cdx2 was increased 5.3- and 12.5-fold in Gsk-3 DKO ESCs 
and S33A ESCs, respectively, but remained unchanged in pi 10* 
ESCs (Figure 3). Finally, qPCR showed that Ido2 expression was 
increased 81.4-fold in Gsk-3 DKO ESCs and 20-fold in S33A ESCs, 
while showing a minimal 1.4-fold increase in expression in pi 10* 
ESCs (Figures). 

Next, we chose to perform similar qPCR validation on genes 
whose expression was changed in Gsk-3 DKO and p 1 1 0* ESCs, but 
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not in S33A ESCs. We selected six genes for analysis - Wnt6, AnxaS 
(Annexin 8), Gata6 (GATA binding protein 6),Aqp8 (Aquaporin 8), 
Foxql (Forkhead box Ql ), and Acta2 (actin, a2, smooth muscle, 
aorta). All of the genes that we assayed showed lower changes in 
gene expression compared to that seen by microarray. For example, 
microarray data showed that Wnt6 and AnxaS were up-regulated 



Table 4 | Summary of comparison between microarray and qPCR gene 
expression data in Gsl<-3 DKO, S33A, and p110* ESCs for genes 
described In this worl<. 





DKO 




S33A 




pllO* 




Microarray 


qPCR 


IVIIcroarray 


qPCR 


Microarray 


qPCR 


Axin2 


6.7x 


7.0x 


3.6x 


3.7x 


1.7x 


-2.5x 


Brachyury 


162. 5x 


622 X 


41. 7x 


38x 


-1.5x 


-1.4x 


Bhmt2 


5.7x 


13.1x 


60. 8x 


56. 5x 


-2.6x 


1.3x 


Cdx2 


27.6 


5.3x 


33.3x 


12. 5x 


1.2x 


1x 


Bhmtl 


6.5x 


14.1x 


49.1 X 


70.1 X 


-5.5x 


2.8x 


ldo2 


9.4x 


81.4x 


3.4x 


20x 


1.3x 


1.4x 


Wnt6 


2.3x 


1.2x 


1.5x 


-1.4x 


10.1 X 


5.7x 


AnxaS 


1.5x 


1.7x 


1.5x 


1.4x 


6.8x 


3.6x 


Gata6 


-1.4x 


-1.1 X 


-3.1 


-2.5x 


-31 .8x 


3.4x 


Aqp8 


-1.5x 


1.6x 


-2.6 


-3.3x 


-48.9 X 


1.7x 


Foxql 


-1.4x 


-l.lx 


-2.9x 


1.5x 


-31 .5x 


-1.7x 


Acta2 


-400.3x 


-5.0x 


-1.8x 


2.3x 


-165.6x 


-2.5x 



Decreases in gene expression are denoted by negative values. 



in pllO* ESCs 10.1- and 6.8-fold, respectively; however, qPCR 
revealed that the expression of Wnt6 and Anxa8 was about half 
of that - 5.7- and 3.6-fold, respectively (Figure 4). Similarly, 
while Gata6, Aqp8, Foxql, and Acta! all showed at least 30-fold 
reductions in gene expression in pllO* ESCs by microarray, the 
qPCR results revealed more subtle changes, ranging from a 3.4- 
fold increase for Gata6 to a 2. 5 -fold reduction in ActaZ (Table 3). 
Importantly, the expression of these genes in S33A ESCs was quite 
different. For example, while Gata6 expression was up 3.4-fold 
in pllO* ESCs, Gata6 expression was decreased 2.5-fold in S33A 
ESCs (Figure 5). Similarly, Acta2 expression was decreased 2.5- 
fold in pllO* ESCs, while Acta2 expression was increased 2.3-fold 
in S33A ESCs (Figure 5). These results are consistent with oppos- 
ing effects of PI3K-mediated insulin signaling and Wnt signaling 
having distinct effects on gene expression. 

DISCUSSION 

The analysis of gene expression changes in Gsk-3 DKO ESCs 
revealed almost 3500 genes whose expression increased or 
decreased at least twofold, reinforcing the notion that Gsk-3 activ- 
ity broadly regulates the expression of many genes. Our data, which 
examined global changes in gene expression in several clonal cell 
lines, illuminated numerous interesting and unexpected features 
of Gsk-3-dependent signal transduction. 

We had initially expected to find many known Wnt target 
genes to be up-regulated in a Gsk-3 and p-catenin-dependent 
fashion and this was indeed true for Axin2 (48) and Brachyury. 
However, most of the other Wnt target genes did not show signif- 
icant increases in expression, and in fact, almost twice as many 



DKO 



S33A 



p110* 





* 






^ 7.0x 
































i 







* 








3.7x 





































9d. 






1 










* 




-2.5x 









Axin2 




FIGURE 2 I Validation of the expression of the Wnt target genes, 
Axin2 and Brachyury In Gsk-3 DKO ESCs, p-catenin S33A ESCs, and 
pllO* ESCs. (A) Quantitative PGR of Axin2 sliowing tlie expression, 
relative to WT control. (B) Quantitative PGR of Brachyury showing the 
expression, relative to WT control. For both results: RQ = relative 




Brachyury 



quantification (all values were normalized to a Gapdti endogenous 
control). Error bars represent standard error of the mean (SEM) between 
biological replicates (n = 3) and technical replicates (n = 3; n = 9 total for 
each gene in each cell line) (*p< 0.001, two-tailed f-test). Fold changes 
compared to WT are shown. 
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FIGURE 3 I Validation of genes whose expression is increased in Gsk-3 
DKO and p-catenin S33A ESCs, but not pllO* ESCs (A) qPCR of Bhmt2 
showing the expression, relative to WT control. (B) qPCR of Cdx2 showing 
the expression, relative toWT control. (C) qPCR of Bhmtl showing the 
expression, relative toWT control. (D) qPCR of ldo2 showing the expression. 



relative to WT control. For all results shown: RQ = relative quantification (all 
values were normalized to a Gapdh endogenous control); error bars represent 
standard error of the mean (SEM) between biological replicates (n = 3) and 
technical replicates (n = 3; n = 9 total for each gene in each cell line) 
{*p < 0.001 , two-tailed f-test); and fold changes compared to WT are shown. 



genes showed decreases in expression. It is possible that this 
simply reflects differences in the regulation of gene expression 
between cell types; alternatively, this data could provide a read- 
out of affected genes due to the constitutive activation of Wnt 
signaling arising from the genetic deletion of Gsk-3a and Gsk-3^. 
Furthermore, since these ESCs have been deficient in Gsk-3a and 
Gsk-3^ since their successful homologous recombination event, 
it is possible that downstream negative feedback effects serve to 
modify the expression of Wnt target genes, making the changes in 
gene expression not appear as dramatic as expected. 

We do, however, identify several putative novel Wnt target 
genes, whose expression is increased in both Gsk-3 DKO and 
P-catenin S33A ESCs. Cdx2, Ido2, Bhmtl, and Bhmt2 were all 



confirmed by qPCR to conform to the microarray expression 
profiling. Interestingly, the transcription factor Cdx2, was recently 
shown to be a direct target of Wnt signaling in the mouse, inde- 
pendently confirming our observation (49). Furthermore, while 
Idol has been shown to be a Wnt target gene (50), Ido2 has not. In 
the mouse genome, the Idol and Ido2 genes lie next to each other 
on chromosome 8 (8qA2) and are transcribed in the same direc- 
tion. Thus, one possibility is that P-catenin-mediated transcription 
could influence both genes from the Idol promoter. Alternatively, 
there maybe additional Lef/Tcf binding sites in the Ido2 promoter. 
An in silico search (51, 52) for Lef/Tcf binding sites in the mouse 
Ido2 promoter identified eight putative binding sites within 2 kb 
of the transcription start site (TSS). Further detailed studies will 
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DKO 



P110* 



S33A 



1.2x 



5.7x 





I 














I 

-1.4x 










L 









Wnt6 



1.7x 



3.6x 




Anxa8 




3.4x 



-2.5x 



Gata6 



1.6x 



1.7x 



* 

-3.3x 



Aqp8 



FIGURE 4 I Validation of genes whose expression is increased in p110* 
ESCs, but not S33A ESCs. (A) qPCR of Wnt6. (B) qPCR of AnxaS. (C) qPCR 
of Wnt6. (D) qPCR of AnxaS. For all results shown: RQ= relative 
quantification (all values were normalized to a Gapdh endogenous control); 



error bars represent standard error of the mean (SEIVI) between biologica 
replicates (n = 3) and technical replicates (n = 3; n = 9 total for each gene in 
each cell line) (*p < 0.001 , two-tailed t test); and fold changes compared to 
WT are shown. 



be required to verify whether any of these putative binding sites 
are bona fide Lef/Tcf binding sites. 

Similarly, Bhmtl and Bhmt2 are also contiguous in the mouse 
genome (13qC3). Neither gene had previously been shown to be 
regulated by Wnt signaling, but we find both genes have increased 
expression in Gsk-3 DKO ESCs (14.1- and 13.1-fold, respectively) 
and even more potently expressed in fi-catenin S33A ESCs (up 
70.1- and 56.5-fold, respectively). Bhmtl and Bhmtl both encode 
for betaine-homocysteine S-methyltransferase enzymes, partici- 
pating in the regulation of homocysteine levels (53). Interestingly, 
variations in BHMTl and BHMT2 have also been implicated in the 
development of cleft palate in certain human populations (54-56). 



Interestingly, mice deficient in GsA:-3f5 develop cleft palate (21, 
24, 25). It is unclear whether this phenotype is a consequence of 
altered expression of Bhmtl or Bhmt2, but based on the data pre- 
sented here, there may be a functional connection between these 
observations. 

While not as dramatic as the changes seen for putative Wnt 
target genes, we found several genes whose expression was insen- 
sitive to constitutively active Wnt signaling, but were increased 
or decreased in expression in a Gsk-3/PI3K-dependent manner. 
It has been demonstrated in animal models of various cancers 
that the transcription factor ZEBl activates PI3K signaling (57- 
59) and this results in increased expression of Gata6 (60). Our 



www.frontiersln.org 



August 2014 | Volume 5 | Article 133 | 9 



Bartman et al. 



Gsk-3-dependent gene expression 



DKO 



D110* 



S33A 











-1.1x 






























1 
















* 




-5.0x 





T 










I 




-1.7x 

















1.5x 








r 





























Foxq1 





















-2.5x 









3 




* 




2.3x 























Acta2 



FIGURE 5 I Validation of genes whose expression is decreased in 
pllO* ESCs, but not S33A ESCs. (A) qPCR of Foxql. (B) qPCR of 
Acta2. For both results: RQ = relative quantification (all values were 
normalized to a Gapdh endogenous control); error bars represent 



standard error of the mean (SEM) between biological replicates (n = 3) 
and technical replicates (n = 3; n = 9 total for each gene in each cell line) 
{*p < 0.001 , two-tailed f test); and fold changes compared to WT are 
shown. 



data from mouse ESCs stably expressing a constitutively active 
form of the pi 10 subunit of PI3K corroborates these findings 
and show that the PI3K-mediated increase in GataS expression 
is not limited to cancer cell populations. Wnt6, AnxaS, and Aqp8 
also showed increased expression in pi 10* ESCs, while all but 
AnxaS had reduced expression in S33A ESCs. While PI3K sig- 
naling has been shown to regulate the subcellular localization of 
Aqp8 in hepatocytes (61), this is the first report that activation of 
PI3K signaling has an effect on Aqp8 transcription. AnxaS expres- 
sion has been shown to be down-regulated by epidermal growth 
factor (EGF) -stimulated PI3K signaling in a model of metastatic 
cholangiocarcinoma (62), while we find that AnxaS expression is 
increased 3.6-fold in pi 10* ESCs, likely reflecting ceU-type differ- 
ences in AnxaS gene regulation. It should be mentioned that we 
selected GataS, Wnt6, AnxaS, and AqpS for validation by qPCR 
because these genes showed increased expression in both Gsk-3 
DKO and pi 10* ESCs by microarray; however, the qPCR results 
showed only modest changes in gene expression in the Gsk-3 DKO 
ESCs, highlighting the importance of not relying on microarray 
gene expression data alone. 

Finally, we also found that the expression of Foxql and Acta2 
were decreased in pi 10* ESCs. Acta2 showed the most robust 
changes in expression (down fivefold in Gsk-3 DKO ESCs, down 
2.5-fold in pi 10* ESCs, and up 2.3-fold in S33A ESCs). We did 
find that Foxql expression is modestly decreased in pi 10* ESCs 
and increased 1.5-fold in S33A ESCs. Foxql has recently been 
shown to be one of the most highly expressed genes in human 
colorectal cancer and has been shown to be a direct target of 
Wnt signaling, which is often constitutively activated in colorec- 
tal cancer cells (63). Based on these data, we expected to find 



increased Foxql expression in both Gsk-3 DKO and S33A ESCs, 
but we did not observe these expected increases. It will be inter- 
esting to determine the mechanism responsible for the lack of 
response by Foxql to constitutively activated Wnt signaling in the 
context of mouse ESCs. 

Taken together, our data shows that, while many genes have 
their expression changed in Gsk-3 DKO ESCs, only a fraction of 
the changes we observe are due to activated Wnt signaling or acti- 
vated PI3K signaling, meaning that there are likely several other 
signaling pathways whose activation affects G5/c-3-dependent gene 
expression. Thus, our data provide a framework for future analyses 
to make these connections. Furthermore, while these studies were 
performed in mouse ESCs, we believe that the findings from this 
study can be used in future studies examining the role of insulin 
or Wnt signaling pathways in other biological settings, such as the 
study of insulin resistance in the diabetic brain. 
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Datasheet 1 | Microarray gene expression data for Gsk-3 DKO ESCs, 
^-catenin S33A ESCs and pllO* ESCs. Tabs at the bottom of the spreadsheet 
contain the microarray data, either as All Data, or sorted by fold changes in 
expression as seen in Gsk-3 DKO, pllO*, or p-catenin S33A ESCs. Within each 
tab is found the Agilent probe ID, Gene symbol. Gene name, NCBI accession 
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number, Entrez Gene ID, fold-changes (FC) compared to WT ESCs, 
false-discovery rates (FDR) compared toWT ESCs, and adjusted P-values (AdjP) 
compared toWT ESCs. The normalized microarray values for each individual 
sample, as well as the averages, are also shown. 
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